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Earlier work has shown that ring dark solitons in two-dimensional Bose-Einstein condensates are 
generically unstable. In this work, we propose a way of stabilizing the ring dark soliton via a radial 
Gaussian external potential. We investigate the existence and stability of the ring dark soliton 
upon variations of the chemical potential and also of the strength of the radial potential. Numerical 
results show that the ring dark soliton can be stabilized in a suitable interval of external potential 
strengths and chemical potentials. We also explore different proposed particle pictures considering 
the ring as a moving particle and find, where appropriate, results in very good qualitative and also 
reasonable quantitative agreement with the numerical findings. 


I. INTRODUCTION 

Over the past few years, there has been an intense 
research interest, not only theoretically, but also exper¬ 
imentally, in the physics of atomic Bose-Einstein con¬ 
densates (BECs) Q, [ 3 , and particularly in the study of 
nonlinear waves Q. Bright M, dark 0 andgap Q 
matter-wave solitons, as well as vortices [3, 0 soli- 
tonic vortices and vortex rings are only some among 
the many structures studied (including more exotic ones 
such as Skyrmions [l^ or Dirac monopoles [l^). 

One of the most prototypical excitations that have 
been intensely studied in experiments are dark soli¬ 
tons 0 . While the early experiments in this theme were 
significantly limited by dynamical instabilities and ther¬ 
mal effects EMI, more recent efforts have been signifi¬ 
cantly more successful in generating and exploring these 
structures. By now, the substantial control of the gener¬ 
ation, and dynamical interactions of such structures has 
led to a wide range of experimental works monitoring 
their evolution in different settings E3“I13|- 

The instability of dark solitons in higher dimensions 
(towards bending [1^ and eventual snaking towards vor¬ 
tices/vortex rings El HI) has been one of the key rea¬ 
sons for the inability to study such states in higher di¬ 
mensions. Although external stabilization mechanisms, 
e.g., utilizing a blue-detuned laser beam [ 2 ^, have been 
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proposed, importantly also variants of such dark solitons 
have been explored in higher dimensions in the form of 
ring dark solitons (RDSs). These efforts were at least 
in part motivated by works in nonlinear optics, where 
they initially were introduced in Ref. [l^, and studied 
in detail, both theoretically (in conservative [28l - l^ and 
—more recently— in dissipative [3l| settings) and exper¬ 
imentally [s^. In turn, RDSs in BEGs were originally 
proposed in Ref. [s^ and their dynamics was analyzed 
by means of the perturbation theory of dark matter-wave 
solitons 0. In other works, RDSs were studied by differ¬ 
ent approaches, e.g., in a radial box [s^, by using a quasi¬ 
particle approach [35| . or by considering them as exact 
solutions in certain versions of the Gross-Pitaevskii equa¬ 
tion (GPE) [ 3 ^. Proposals for the creation of RDS, e.g., 
by means of BEG self-interference or by employing 
the phase-imprinting method [s^ . as well as generaliza¬ 
tions of such radial states (including multi-nodal ones) 
[a [11 have also been considered. Moreover, generaliza¬ 
tions of RDSs were studied in multi-component settings, 
in the form of dark-bright ring solitons [i^l (emulating 
the intensely studied context of multi-component one¬ 
dimensional (ID) dark-bright solitons pll - l44l |). or in the 
form of vector RDS in spinor E = 1 BECs [11] • 
portantly, structures of the form of radially symmetric 
dark solitons, closely connected to RDSs exist also in 
three-dimensions with a spherical rather than cylindrical 
symmetry (so-called “spherical shell solitons” 

Nevertheless, in none of these contexts (either one- or 
multi-component), was it possible to achieve complete 
stabilization of the RDSs. In particular, stabilization 
mechanisms that have been proposed, e.g ., by “filling” 
the RDS by a bright soliton component |4C| or by employ¬ 
ing the nonlinearity, management (alias “Feshbach reso¬ 
nance management” [^ ) technique [d^, were only able 
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to prolong the RDSs’ life time. In fact, it was illustrated 
that the instabilities of the ring-shaped solitons were 
connected, bifurcation-wise, to the existence of vortex 
“multi-poles”, such as vortex squares (which are generi- 
cally stable in evolutionary dynamics), vortex hexagons, 
octagons, decagons etc.; all of these states are progres¬ 
sively more unstable. This picture has been corroborated 
by detailed numerical computations in Ref. . It is our 
aim in this work to revisit the RDSs and their destabiliza¬ 
tion mechanisms and, indeed, to propose a technique for 
their complete dynamical stabilization. Our technique is 
reminiscent of that of Ref. [1^ in that we introduce a 
potential induced by a radial blue-detuned laser beam. 
Radial potentials of a similar form have been intensely 
used in recent experiments, e.g., by the groups of [48| 
and [ 4 ^ and are hence accessible to state-of-the-art ex¬ 
perimental settings. 

Our presentation of this effort to stabilize the RDS 
in the form of a dynamically robust state of quasi-two- 
dimensional BECs can be summarized as follows: we 
introduce, in Sec. im the mathematical model and our 
specihc proposal towards a potential stabilization of the 
RDS. We also incorporate in this Section theoretical at¬ 
tempts to explore the coherent dynamics of the ring soli- 
ton by means of a particle model. Our numerical results 
are presented in Sec. nni initially revisiting (for reasons 
of completeness and to facilitate the exposition) the case 
without the external radial barrier potential and subse¬ 
quently incorporating it in the picture. Finally, our con¬ 
cluding remarks are presented in Sec. IIVI and a number 
of important open future directions is also highlighted. 


II. MODEL AND MATHEMATICAL SET-UP 

A. The Gross-Pitaevskii equation 

In the framework of lowest-order mean-field theory, 
and for sufficiently low-temperatures, the dynamics of 
a quasi-2D (pancake-shaped) BEC confined in a time- 
independent trap V (r) is described by the following di¬ 
mensionless GPE 

+ IV'PV'- (1) 

where '0(a:, y,t) is the macroscopic wavefunction of the 
BEC, fi is the chemical potential, and V{r) (with r = 
\/+ y^) is the external potential. The latter, is as¬ 
sumed to be a combination of a standard parabolic (e.g., 
magnetic) trap, I4it(?'), and a localized radial “pertur¬ 
bation potential”, Upert (r ), namely: 

V{r) = Vmt(?') + Vpert(?') = + ^pert(T’), (2) 

with D being the effective strength of the magnetic trap. 
For the numerical results in this work we chose a nominal 


value of D = 1 unless stated otherwise. As will be evi¬ 
dent from the scaling of our findings below, the particular 
value of D will not play a crucial role in our conclusions. 

The GPE in the Thomas-Fermi (TF) limit of large /i 
has a well known ground state '0 tf = \/ max(/i — U, 0). 
The other interesting limit is the linear one where the 
self-interaction term can effectively be ignored. In this 
limit, the GPE reduces to the 2D harmonic oscillator 
problem. Both limits are particularly useful for our con¬ 
siderations: the former enables the consideration of the 
ring-shaped soliton as an effective particle, the latter en¬ 
ables the construction of the ring as an exact solution 
in the linear limit, which is continued in the nonlinear 
regime. 

Here, we will focus on the single RDS which, in the 
linear limit, can be viewed as a superposition of the |2,0) 
and |0, 2) quantum harmonic oscillator states, namely: 

l; \ |2, 0) -|- |0, 2) 2 -nr^/2 ,'o^ 

|V'RDS)linear = -- OC (r - l)e ' . (3) 

v2 

This linear state, which exists for /i > SD (i.e., beyond 
the corresponding linear limit of the above degenerate 
n-|-m = 2 states, where n and m are the respective indices 
along the x- and y-directions, characterizing the quan¬ 
tum harmonic oscillator sate |n, m)), can be continued to 
higher chemical potentials. However, the RDS is known 
to be inherently unstable for all values of fi beyond the 
linear limit [13,113113 ■ This instability breaks the orig¬ 
inal radially symmetric state into vortex multi-poles, as 
originally shown in Ref. [s^ and subsequently examined 
from a bifurcation perspective in Ref. [^. Our scope is 
to provide a systematic understanding of the RDS insta¬ 
bility modes and how to suppress them, so as to poten¬ 
tially enable its experimental realization. Similar consid¬ 
erations in the context of exciton-polariton condensates 
(where a larger range of tunable parameters exists due to 
the open nature of the system and the presence of gain 
and loss) have led both to the theoretical analysis [5l| 
and to the experimental observation of stable RDSs. 

Following the motivation of the earlier work of Ref. 
on planar dark solitons, in conjunction with the recent ex¬ 
perimental developments in the context of radial [i^ and 
more broadly, in principle arbitrary, so-called painted [49j | 
potentials, we propose the following form for Ipert('c): 

(4) 

where Tc, A and a represent, respectively, the radius, the 
amplitude and the width of this ring-shaped potential. 

Since RDSs feature radial symmetry, we first express 
Eq. ([1} in the form: 

+ Vir)^|; + (5) 

We also assume that a stationary RDS state, ip = ip{r, t), 
governed by the effectively ID model (O, is characterized 
by a radius r^. In other words, we will hereafter opt to 
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locate the perturbation potential at the fixed equilibrium 
position of the RDS. For our analysis, the control param¬ 
eters will be the strength A of the perturbation potential 
and the nonlinearity strength (characterized by the chem¬ 
ical potential /r); as concerns the width tr of V^ert(?'), it 
will be fixed (unless otherwise stated) to the value cr = 1, 
which is of the order of the soliton width —i.e., of the 
healing length. 

Below, we proceed with the study of the effect of the 
perturbation potential on the existence and stability of 
the RDS. Stability will be studied from both the spec¬ 
tral perspective, through a Bogolyubov-de Gennes (BdG) 
analysis, and from a dynamical time evolution perspec¬ 
tive. The latter, will involve direct numerical integra¬ 
tion of Eq. o, whereby a (potentially perturbed) RDS 
is initialized and its evolution is monitored at later times. 
On the other hand, BdG analysis for a stationary RDS, 
i/joif)-, will involve the study of the eigenvalue problem 
stemming from the linearization of Eq. o, upon using 
the perturbation ansatz: 


the system, take the following form 
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The above system suggests the existence of stationary 
RDSs, due to the interplay (to the leading-order approx¬ 
imation in D) of an effective attractive trapping potential 
and an effective curvature-induced repulsive logarithmic 
potential —see first and second terms in the right-hand 
side of Eq. ®, respectively. A more systematic analy¬ 
sis, that takes into regard higher-order terms in D, shows 
that the critical radius for which a stationary ring exists 
is given by 


ipix, y, t) = ipoir) + S (^u{x, y)e^* -h v*{x, ?/)e^** j , (6) 


where [A, (u, w)"^] is the eigenvalue-eigenvector pair, S is 
a formal small parameter, and the asterisk denotes com¬ 
plex conjugation. Then, the existence of eigenvalues with 
non-vanishing real part signals the presence of dynamical 
instabilities. These come in two possible forms: (a) gen¬ 
uinely real eigenvalue pairs, which are associated with an 
exponential instability; and (b) complex eigenvalue quar¬ 
tets that denote an oscillatory instability, where growth 
is coupled with oscillation. The above symmetry of the 
eigenvalue pairs (i.e., the fact that they only arise in pairs 
or quartets) stems from the Hamiltonian nature of the 
problem. 


B. The particle picture for the ring dark soliton 

A natural way to obtain a reduced dynamical descrip¬ 
tion of the RDS is to adopt a particle picture and use a 
variational approximation discussed in detail in Ref. . 
According to this approach, in the TF limit (i.e., for 
sufficiently large chemical potential), the RDS state can 
be approx imated by a pr oduct of the TF ground state, 
V’TF = Y^max(/r — E, 0), and a (potentially traveling) 
dark soliton of radial symmetry, of the form: 

ipDsirA) = b{t)ta.nh.[y/JIb{t)(r - rc{t))] +ia{t), (7) 

where b and a (with + b'^ = 1) set, respectively, the 
depth and velocity of the soliton, while Vc is the RDS 
radius. Then, the Euler-Lagrange equations for the two 
independent effective variational parameters Xc and a, 
stemming from the averaged renormalized Lagrangian of 


r. = (10) 

Notice that, according to the discussion of Ref. and 
in accordance with the computational analysis presented 
below (see Sec. lIIIG]) . the numerical results strongly sug¬ 
gest an asymptotic critical radius Xc = -\/^/2/D (see also 
the discussion in Refs. [H, HI, [b^). 

This discrepancy suggests the consideration of alter¬ 
native ways of determining the stationary RDS’ radius. 
Here, for reasons of completeness, we will present such 
an alternative approach, based on the earlier work of 
Ref. [ij for a different system (namely, ring-like steady 
state solutions of coupled reaction-diffusion equations). 
More specifically, our starting point will be the steady 
state problem associated with Eq. ([5]), where we will 
“lump” the potential terms as V{x) = HMT+Epert(?')- Us¬ 
ing the ansatz '4’{x) = 'ipTpir)q{x), we obtain the steady 
state problem: 

y' + ^,q{l-q^) = P{r), ( 11 ) 


where 

P{x) = Vqil -<Y)-Yr 


'4’tF V'tF / 1 V'tf 


and primes denote derivatives with respect to x. Then, 
seeking a stationary RDS solution in the form of q{x) = 
t&vth{y/Ji{x — Xc)) and multiplying both sides by q' in 
Eq. (fTl]) . we find that the left-hand side is simply dH/dx^ 
where PI is the effective Hamiltonian PI = g'^/4 — /x(l — 
q^Y!A. Hence, upon integrating in x from —oo to oo, 
bearing in mind that the error between r = 0 and 
r —> — oo is exponentially small, we obtain the explicit 
solvability, Melnikov-type, condition [b^ : 


/ OO 

P{x)q {x)dx = 0. 

-OO 


( 12 ) 
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Upon evaluating the integrals of all five terms associated 
with P{r) within Eq. (TT^ . we should obtain an algebraic 
equation for the equilibrium position of the RDS. Indeed, 
evaluating the first potential term (for ^ = 0), through 
a series of rescalings and integrations by parts, leads to 
V.'^rc/{3y^). In turn, the second term yields —2^j (3rc) 
and the fourth term yields 2fl'^rc^/JI/[3{fJ, — f2^r^/2)], 
while the other terms contribute at higher order. Putting 
all the terms together in the case of A = 0 yields the pre¬ 
diction 


rc = 



(13) 


where a = 4 — 2\/3 ~ 0.5359; this result is more accurate 
than the one of Eq. m. as will be discussed in more 
detail in Sec. IIII Cl 

Finally, we proceed to give a third method, based on 
the analysis of Ref. [s^, that will prove to be the most 
accurate one in connection to our computations of not 
only statics but also dynamics of RDS states in the nu¬ 
merical section that will follow. In the latter approach, 
it is argued that the equation of motion can be derived 
by a local conservation law (i.e., an adiabatic invariant) 
in the form of the energy of a dark soliton under the 
effect of curvature and of the density variation associ¬ 
ated with it. More specifically, knowing that the en¬ 
ergy of the one-dimensional dark soliton is given by Q 
^DS = (4/3)(/i - aic)^/^, where Xc is the dark soliton 
position, the generalization of the relevant quantity in 
a two-dimensional domain bearing density modulations 
reads: 



FIG. 1: (Color online) Top panels: the RDS’ real-valued 
profile (left) and the corresponding density plot (right) for 
= 16. Middle left: a radial profile of the relevant state. 
Middle right: number of particles N = J dr as a function 
of chemical potential fj,, showing the continuation of states 
from the linear limit to the nonlinear regime. Bottom pan¬ 
els: the imaginary (left) and real parts (right) of the spec¬ 
trum; showcased is the generic instability of the RDS, and 
the emergence of additional unstable eigenmodes thereof as 
is increased. The value for the trap strength in this figure 
and all remaining figures is D = 1 (unless stated otherwise). 


III. RESULTS 


fRDs(?') = 27r7’ 


^(/r-U(r)-f2)3/2 


(14) 


Thus, by assuming this quantity is constant, namely 
^RDs(^) = ^RDs(u), where Tc is the equilibrium loca¬ 
tion of the ring, we obtain an equation for f^. Taking 
another time derivative on both sides, we finally obtain 
Newtonian particle dynamics for the ring in the form: 


r 


1 W 

2 ~^ 


1 

3r 





(15) 


When ^ = 0, this equation of motion for the RDS posi¬ 
tion yields the equilibrium = y/ /j,/2/D, a result which, 
as highlighted also above and as will be demonstrated 
below, is the one most consistent with the numerical ob¬ 
servations. This, in turn, motivates us to use the above 
approach of Ref. [s^ not only for the statics, but also for 
the dynamics in the following section and additionally, 
not only for the case without the radial defect of A = 0, 
but also for that bearing the radial defect i.e., for A ^ 0. 

We now proceed to test these predictions, as well as 
to examine the BdG stability analysis and the dynamical 
evolution of the RDS, both in the absence (initially, for 
comparison and guidance) and then in the presence of 
the radial perturbation potential. 


First, we briefly summarize the numerical techniques 
used in this work. Stationary states in both ID (i.e., in 
a radial form) and 2D were identified using a centered 
finite-difference scheme within Newton’s method. The 
spectrum of the stationary states (i.e., the result of the 
BdG analysis) was calculated using the eigenvalue prob¬ 
lem derived from Eq. Finally, for the dynamics of 
the system, we used direct integration employing second 
order finite differences in space and fourth-order Runge- 
Kutta in time. 


A. Basic properties of the ring dark soliton 

Let us start by summarizing some of the basic prop¬ 
erties of the RDS without the perturbation potential. A 
typical RDS state in the TF limit of large chemical poten¬ 
tial /i is shown in the top and middle left panels of Fig. [T] 
the top right panel shows the corresponding density. As 
indicated in the previous section, the RDS has a linear 
limit (built out of the eigenstates of the 2D quantum 
harmonic oscillator). The continuation of such a state in 
the nonlinear regime is shown in the middle right panel 
of Fig. [1] The imaginary and real parts of the spectrum 
of the RDS are shown in the bottom panels of Fig. [T] 
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FIG. 2: (Color online) The most unstable modes at a few 
representative chemical potential values — 4 , 6 , 9, 11, 14, 
and 16 (from left to right, top to bottom) associated with the 
instability of the RDS. Left and right subpanels correspond, 
respectively, to the absolute value and phase of the modes. 
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FIG. 3: (Color online) Dynamics ensuing from the unstable 
RDS for — 16. Note that the RDS first deforms into seven 
pairs of vortices (in accordance with the most unstable mode 
for these parameters values; see the bottom right panel of 
Fig.E}, and then eventually turns into a dynamical evolving 
vortex cluster for longer times. During evolution, some of the 
vortices are “absorbed” by the BEG periphery and the system 
is eventually left with four interacting vortices. 


Note that the RDS is unstable for any value of /r beyond 
the linear limit. More importantly, in line with what was 
also presented in Ref. [^, as /i increases, more unstable 
modes keep emerging, through eigenvalue pairs that cross 
through the origin. These signal pitchfork bifurcations, 
to which we now turn. 

Studies of RDS in atomic BECs have illustrated their 
dynamical breakup into vortex-antivortex pairs (see, e.g.. 
Refs. [11,13^). To complement this picture, we now dis¬ 


cuss the most unstable modes of the BdG analysis. Some 
representative eigenmodes at /i = 4, 6, 9, 11, 14 and 16 
are shown in Fig. [2] It is interesting to observe that the 
identified modes indicate a clear connection to an increas¬ 
ing number of pairs of vortices. The first unstable mode 
appears to be connected to two-pairs, i.e., to a vortex 
quadrupole. Indeed, the vortex quadrupole exists as a 
state |56l | for any value of fx beyond the linear limit of 
/i = 3D, being constructed as: 


I linear 


|2,0)+f|0,2) 

72 


(16) 


Subsequent destabilization modes reveal a three¬ 
fold symmetry (leading to the bifurcation of vortex 
hexagons 0), a four-fold symmetry (leading to vor¬ 
tex octagons), then a five-fold (decagons), a six-fold (do¬ 
decagons), and so on. These different eigenvectors are 
clearly illustrated in Fig. [2] and the existence and sta¬ 
bility of the corresponding emerging (from the pitchfork 
bifurcation) vortex n-gon cluster states is discussed in 
Ref. [ 53 . 

A dynamical study of the states shows that the evo¬ 
lution initially results in vortex pairs, in agreement with 
Fig. [2] However, gradually some vortices may move out 
of the BEG and get lost in the background, leaving be¬ 
hind a complex, interacting cloud of vortices, as shown 
for ^ = 16 in Fig. The resulting interaction dynamics 
between vortices in the cluster, and the associated trans¬ 
fer of energies between different scales, may represent 
a very interesting setting for exploring turbulence phe¬ 
nomena and associated cascades in line, e.g., with recent 
experimental efforts of Ref. . 


B. Adding the perturbation potential 

Having analyzed the unperturbed case, we now exam¬ 
ine the case with the radial Gaussian potential. The ex¬ 
istence of the RDS structure in the latter case can be 
captured as a function of (A,/i) —see Fig. |4l We used 
maxdTl) (i.e., the max root density) as a diagnostic in¬ 
stead of N for practical visualization purposes, in this 
case. We can see that for a fixed value of the den¬ 
sity decreases as A increases (a natural feature, given 
the repulsive nature of the perturbation potential) until 
a critical value of A —shown as a purple line— is reached, 
beyond which the RDS will cease to exist. In the linear 
limit of ^ = 3D, even a very small positive perturbation 
of A will destroy the RDS state. The monotonic depen¬ 
dence of on the critical A appears to be approximately 
linear. 

We proceed now with the central theme of this study, 
which is the dynamical stabilization of the RDS. To char¬ 
acterize the stability of the RDS in the {A, /r) plane, in 
Fig. [S]we show a plot of the max(Re(A)) as a function 
of {A,fi). The right most purple line, as before, depicts 
the critical values of A beyond which no RDS solution 
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FIG. 4: (Color online) Max(|4'|) as a function of {A,fj,). Note 
that N decreases as A increases when holding fixed until 
some critical set of values of A (depicted by the purple line) 
beyond which the RDS will cease to exist. In the linear limit 
fj, — 3II, even a very small positive perturbation of A will 
destroy the RDS state. 



FIG. 5: (Color online) Instability growth rate max(Re(A)) 
as a function of {A, /r). The area between the left (green) 
and right (purple) curves corresponds to the region where the 
RDS exists and with vanishing max(Re(A)), i.e. the RDS 
is completely stable. The rightmost purple line is also the 
boundary of the critical values of A beyond which no RDS 
solution exists. 


exists. The region enclosed between the green and pur¬ 
ple lines corresponds to the regimes where RDS exists 
with vanishing max(Re(\)), i.e., the RDS is completely 
stabilized by the presence of the external Gaussian ring 
perturbation potential. One interesting feature is that 
the relevant stability landscape is rather complex with 



A 

FIG. 6: (Color online) Cross section of the instability growth 
rate max(Re(A)) at /r = 4. The right most point of the curve 
corresponds to the critical value of A beyond which no RDS 
solution exists. The two blue squares are two points in two 
different instability regimes but with similar instability rates 
whose full spectrum is shown in Fig. 0 


potential sequences of destabilization and restabilization 
for values of /i > 3.6 (we will return to this point below). 
However, the principal conclusion obtained from Fig.[S]is 
that the RDS is generically subject to full dynamical sta¬ 
bilization for any value of the chemical potential and for 
suitable intervals of the perturbation potential strength 
A in the vicinity of the linear limit. The feature that 
the stabilization is enabled near the linear limit is rather 
natural to expect also on the basis of our earlier results 
for H = 0 in Fig. [T] Given that the RDS is progressively 
more and more unstable (with a higher number of desta¬ 
bilizing modes) as /r increases suggests that the perturba¬ 
tion potential may be unable to suppress this multitude 
of unstable modes, especially far from the linear limit. 

To gain further insight on this stability plane, let us 
now study a typical cross section of Fig. [5] at /r = 4. The 
cross section is shown in Fig. |6l A detailed study of the 
full spectrum shows the existence of two intervals of in¬ 
stability which are not of the same nature. The leftmost 
interval (including A = 0 in the absence of a defect) cor¬ 
responds to a typically large(r) growth rate. Here, the 
instability derives from real eigenvalue pairs. Gonnecting 
with Fig. [Hand the case of A = 0, we recognize that this 
unstable mode is associated with the breakup to vortex 
quadrupoles. As A becomes increasingly more negative 
to the left of the figure, other modes may, in turn, domi¬ 
nate the instability dynamics (the “bend” in the stability 
diagram represents such a “take-over” of the dominant 
instability by a different mode; cf. Fig. [T|). However, it 
is observed that as A increases on the positive side, the 
unstable real pair(s) decrease in their real part and even¬ 
tually cross through the origin of the spectral plane, be- 
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FIG. 7: (Color online) Full stability spectrum correspond¬ 
ing to the two blue squares in the two different instability 
regimes in Fig. [5] for /r = 4. Left and right panels corre¬ 
spond to j 4 = 1.07 and A — 1.28, respectively. Note that the 
two regimes do not share the same nature of instability. The 
large amplitude case (left) has the instability on the real axis 
(i.e., exponential instability) while the small amplitude case 
(right) has the instability in the form of a complex quartet 
(oscillatory instability). 


coming imaginary and hence stabilizing the RDS state. 
This is, once again, a key finding of our work, repre¬ 
senting the RDS stabilization. However, as the (formerly 
unstable) eigenvalues bear a so-called ‘negative energy’, 
upon climbing up the imaginary axis, they may collide 
with eigenvalues associated with ‘positive ener^’ modes 
(see, e.g., the discussion in pp. 56-58 of Ref. 0). This 
type of collision gives rise to a complex eigenvalue quartet 
and a different (weak) oscillatory dynamical instability, 
or a so-called Hamiltonian Hopf bifurcation; see, e.g., the 
discussion of Ref. . The latter scenario leads to small 
instability bubbles, as the quartet may form, but sub¬ 
sequently the eigenvalues may return to the imaginary 
axis, splitting anew into two imaginary pairs. 

The two (exponential and oscillatory) instability sce¬ 
narios are illustrated in the two panels of Fig. [7] for 
smaller and larger values of A, respectively. The most 
unstable mode of each state is shown in Fig. [51 illustrat¬ 
ing the distinct nature of the instability in the different 
scenarios. The state at A = 1.07 is in the same branch of 
A = —1,0 and 1, and its instability leads to a deforma¬ 
tion towards a vortex quadrupole state in a way similar 
as the hrst plot of Fig. [21 On the other hand, the state 
at A = 1.28 appears to have a different type of instabil¬ 
ity that instead resembles a vibrational mode (the type 
of mode that could be captured through a ring particle 
model). The time dynamics of the two states are shown 
in Fig. [5] and Fig. [TUI respectively. In the former case, we 


FIG. 8: (Color online) The most unstable modes of Fig. [T] 
Top and bottom row of panels correspond, respectively, to the 
absolute value (left subpanels) and phase (right subpanels) of 
the solutions for the left and right cases depicted in Fig. [T] 



FIG. 9: (Golor online) Dynamics of the state in the top panels 
of Fig. (8] The odd panels depict the absolute value of the field 
while the even panels depict its phase. The state is oscillating 
between the vortex quadrupole and the RDS, but very weakly. 


observe the recurrent formation of a vortex quadrupole 
(this is not immediately discernible in the density but dis¬ 
tinctly visible in the phase pattern), in accordance with 
the identified unstable mode. In the latter, indeed un¬ 
stable vibrational dynamical characteristics can be seen 
in the motion of the ring, which, however, appears to 
maintain its radial structure. 

A different cross section of the stability plane of Fig. [5] 
is given in Fig. [TTl now for the case of A = 0.5, and 
varying the chemical potential /i. From this perspective, 
we observe that A delays the onset of instabilities as ^ is 
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FIG. 10: (Color online) The same plot as Fig. [5] but for the 
state in the bottom panels of Fig. [8] This state has a different 
nature of instability from the one in Fig. The instability is 
like a vibrational mode. 



FIG. 11: (Color online) Cross section of max(Re(A)) at A = 
0.5. The solution starts to exist around /r = 3.35. Note that 
A delays the set in of instabilities as is increased. 


FIG. 12: (Color online) Time evolution of states at ^ = 4 with 
A = —1, 0 and 1 (top to bottom rows of panels). Note that 
the state of A = 1 is significantly less unstable than those of 
A = — 1 and A = 0, which have roughly the same instability 
growth rate. Note also that all three states deform toward 
the vortex quadrupole state initially, although the third one 
maintains an oscillatory pattern between a recurring ring and 
a vortex quadrupole. 
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FIG. 13: (Color online) Time evolution of states at ^ = 4 with 
A — 1.14 which is in a completely stable parametric interval. 
The state is shown to be stable up to t=1000, in agreement 
with our spectral results. 


increased. Another way to look at the effects of A and 

is that A plays effectively the opposite role to that of 
/i: the increase of A (for fixed /r) drives the eigenmodes 
away from the real axis and into the imaginary axis while 
the increase of the chemical potential for fixed A drives 
the eigenmodes away from the imaginary axis and into 
the real axis, causing instability. We believe that this 
discussion provides a unified perspective on the sources 
of destabilization and the potential for re-stabilization of 
the RDS. 

In all the cases considered, the stability conclusions 
were also found to be consonant with the corresponding 
dynamics, of which we now present a few additional case 


examples. In particular, we study the dynamical evolu¬ 
tion of states at /i = 4 for different values of A = — 1, 0,1 
(see Fig. [12]) and 1.14 (see Fig. [T3|) to probe the effects 
of the variation of A. Note that the cases of A = — 1 
and 0 are about equally unstable at /i = 4 with A = — 1 
bearing a slightly larger growth rate. The case of A = 1, 
however, is very close to, albeit not within the stabiliza¬ 
tion regime. On the other hand, the case of A = 1.14 is 
fully stabilized. We add a random perturbation to the 
states, ensuring that the number of atoms in each case 
is, upon perturbation, 1.0013 times of the unperturbed 
one. The results of the dynamical evolution of the for- 
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FIG. 14: (Color online) The location of the RDS scaled by 
as a function of (thick solid blue line). Note that the 
numerical values reach a limiting value of ll\/2 (thin horizon¬ 
tal solid red line) when /r is large. The particle picture can 
approximately describe the behavior and over estimates 
Tc , but nevertheless is still an interesting approximate descrip¬ 
tion of the RDS. The particle approach using the perturbed 
Lagrangian method [see Eqs. ([8]) and Q] corresponds to the 
thin dotted-dash green line while the solvability condition for 
the steady state problem method [see Eq. (1 1311 ] is depicted by 
the thin dashed black line. 


mer three cases are shown in Fig. [T^l Note that both 
states for A = —1 and A = Q are relatively quickly de¬ 
formed around t = 25 while the state for A = 1 deforms 
only much later around t = 70, due to its weaker growth 
rate. In all three cases, the states evolve initially towards 
the vortex quadrupole waveform. While the former two 
states will quickly deform afterwards and lose their radial 
structure, the third state can oscillate between the RDS 
state and the vortex quadrupole state for a much longer 
time at least up to t = 1000. A dynamical evolution of 
states at A = 1.14, which is in a completely stable para¬ 
metric interval, is shown in Fig. 1131 The state is shown 
to be stable at least up to t = 1000, in agreement with 
the spectral findings and corroborating the full stabiliza¬ 
tion achieved by the presence of the Gaussian repulsive 
impurity. 


C. The particle picture of the ring dark soliton 

We first study how the equilibrium location Vc of the 
RDS changes with chemical potential /r, especially in the 
large density limit without the perturbation potential, 
and compare the numerical results and the particle pic¬ 
ture predictions. Numerical results (for D = 1) suggest 
that Tc = y/ /i/2 (see thin horizontal red line) in the large 
/j, limit as shown in Fig. [141 As mentioned in Sec. IIIBI 
a systematic analysis of Eqs. ([Sj) and (jH)) yields the es¬ 
timate Tc = y^O.SOIO^/D (see horizontal thin dotted- 
dash green line). On the other hand, using the solv¬ 



FIG. 15: (Golor online) Dark soliton width w as a, function of 
the chemical potential /r. The blue solid line corresponds to 
fitting a profile ipTFir) x tanh(y/w(r — Tc)) to the PDE steady 
state for fl = 0.2. The red dashed line corresponds to the 
approximate value of the background at the location of the 
RDS, see Eq. (|T7||. 


ability condition for the steady state problem described 
in Sec. Ill B[ one obtains the better prediction of the RDS 
position Tc = ^/0.5359fJ,/^2; see Eq. ([T^ and thin hor¬ 
izontal dashed black line in Eig. [TH It is important to 
mention that, although the above two particle approaches 
are able to capture the y/Jt/fl behavior of Tc, they do not 
lead to the precise numerical prefactor. This may be at¬ 
tributed to the choice of the ansatz 0, where the width 
of the stationary dark soliton is chosen to be y/Jl. This 
selection corresponds to the width of a dark soliton in a 
homogeneous background of density fi. However, due to 
the non-homogeneity of the BEG background, the RDS 
placed at Tq experiences a background density which 
can be approximated using the TF regime (valid for large 
/i) to be 

^lo-1|^TF{rc) = ^^-V{rc)=^i-^^'^rl. (17) 

For instance, in Fig. [T5| we show an example where we 
extracted the width of the dark soliton for D = 0.2 as 
a function of /r. As it is clear from the figure, as /r in¬ 
creases, the width of the dark soliton converges to y/]Io 
as prescribed in Eq. (HZl). Lastly, it is relevant to point 
out that, remarkably, the adiabatic invariant theory of 
Ref. [ 3 ^ properly captures the asymptotic growth of the 
radius of the RDS as increases. It is for that reason that 
we will hereafter utilize the particle picture of Eq. m 
and Ref. for our further static and dynamics consid¬ 
erations. 

We now study the effect of A on Tc. Eigure [TBl depicts 
Tc as a function of A for /r = 24. It is clear that the par¬ 
ticle picture can capture the effect of A fairly accurately. 
It is also observed that the critical radius decreases in 
comparison to the A = 0 limit, in the presence of a re- 
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FIG. 16: Position of the RDS as a function of A for ^ = 24. 
The (red) circles correspond to the full PDE dynamics and 
the (green) triangles to the particle picture (PP) described in 
Sec.lITB] 


pulsive defect, while the opposite is true in the case of an 
attractive defect. 

Finally, we study the radial oscillatory motion of the 
RDS in both the case bearing and in that without the 
perturbation potential. We initialize our displaced RDS 
state by superposing a suitable hyperbolic tangent profile 
to (i.e., multiplying it with) the numerically exact ground 
state at the same chemical potential /i = 24. Note that 
the RDS is unstable at such a high chemical potential, 
therefore, we can only simulate the PDE dynamics for a 
limited amount of time, before an instability leading to a 
polygonal cluster of vortices ensues. The comparison of 
the PDE and the particle picture dynamics for the cases 
of A = 0 and A = 1 are shown, respectively, in the top 
and bottom panels of Fig. [T3 We see that the particle 
picture is able to capture the essential PDE radial os¬ 
cillation dynamics both with and without the Gaussian 
barrier. 


IV. CONCLUSIONS AND FUTURE 
CHALLENGES 

In this work, we studied the existence and stability of 
ring dark soliton states, initially in the absence and sub¬ 
sequently in the presence of a radially localized Gaussian 
perturbation potential. We have systematically shown, 
via a combination of spectral analysis and direct numer¬ 
ical simulations, that the ring dark soliton can be stabi¬ 
lized by adding the perturbation potential with a suitable 
strength, for all values of the chemical potential that we 
have considered herein. Our systematic spectral analysis 
has also revealed why this stabilization mechanism can 
only be effective near the linear limit of the system. It 




FIG. 17: Radial oscillatory motion of the RDS with fj, — 24 for 
A = 0 and A = 1. The central radius of the RDS is extracted 
from the PDE dynamics (green dots) and compared to the 
ODE evolution of the particle picture (PP, red line) according 
to Eq. (fTSl) . 


has also revealed the potential for secondary instabilities 
(due to pair collisions on the imaginary axis and com¬ 
plex eigenvalue quartets emerging from them) due to the 
excited state nature of the ring. 

An additional effort, significantly motivated by the po¬ 
tential of the above method to lead to stable RDS vibra¬ 
tions, was that of deriving dynamical equations for their 
motion. We evaluated different techniques to this effect, 
showcasing the fact that although all approaches gave 
fairly similar results, the adiabatic invariant method of 
Ref. [S^ presented a distinct advantage in capturing the 
radius of stationary rings. A self-consistent perturbative 
technique (based on earlier work in reaction-diffusion sys¬ 
tems) was also adopted to that effect and was shown to 
give reasonably accurate results in its comparison with 
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the full numerical results. Going beyond the “station¬ 
ary particle” approach, allowing motion along the radial 
direction, an intriguing goal for the future may be to ex¬ 
amine the ring soliton as a filamentary pattern embedded 
in 2D, which, in addition to radial internal modes, may 
possess bending ones (but without breaking). Such stud¬ 
ies may in turn enable the observation of collisions and 
deformations of rings upon interactions, a topic that has 
been of interest also in nonlinear optics [s^. 

Finally, it may well be relevant to explore settings 
beyond the realm of two spatial dimensions, extending 
the present considerations to the case of 3D solitonic or 
vortex rings and other such patterns. Earlier work es¬ 
tablished how to construct such states in isotropic and 
anisotropic 3D limits starting from linear eigenstates . 
It is then of particular interest to continue such states in 
the nonlinear realm and explore their spectral and dy¬ 
namical stability using tools similar to the ones proposed 
herein. Efforts along these directions are currently in 
progress and will be presented in future publications. 
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